function sc = subclusters(vox, val, sep)
% sepclusters  - separate mega cluster into sub-clusters
%
% FORMAT:       sc = subclusters(vox [, val], sep)
%
% Input fields:
%
%       vox         either a 3D volume with values or a Vx3 list
%       val         required if vox is a list of coordinates
%       sep         cluster separation (in voxel resolution)
%
% Output fields:
%
%       sc          1xC sub clusters structure
%        .peak      peak coordinate
%        .val       values
%        .vox       voxel coordinates

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end

if numel(val) == 1
    sep = val;
    msk = double(vox > 0);
    [vox, voy, voz] = ind2sub(size(msk), find(msk));
    vox = [vox, voy, voz];
else
    vox = round(vox);
    vox = 1 + vox - ones(size(vox, 1), 1) * min(vox);
    msk = zeros(max(vox));
    msk(sub2ind(max(vox), vox(:, 1), vox(:, 2), vox(:, 3))) = 1;
end

% sort values
[sv, svi] = sort(val);

% start with marking highest value as first cluster
cc = 2;
msk(vox(svi(1), 1), vox(svi(1), 2), vox(svi(1), 3)) = cc;

% iterate until all values are used
